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EFFECT OF BOTTOM TOPOGRAPHY, EDDY DIFFUSIVITY, AND WIND 
VARIATION ON CIRCULATION IN A TWO-LAYER STRATIFIED LAKE 
by Richard T. Gedney, Wilbert Lick,* and Frank B. Molls 
Lewis Research Center 

SUMAAARY 

The steady-state, wind-driven circulation is calculated in a stratified lake com- 
posed of two layers having uniform but unequal densities and eddy diffusivities. The 
position of the thermocline and the three dimensional velocities in both layers are cal- 
culated from an asymptotic solution of the shallow lake equations when the Ekman num- 
ber in the epilimnion (upper layer) is of order one but the ratio of hypolimnion (lower 
layer) to epilimnion eddy diffusivities is much less than one. This analysis differs from 
previous ones in that the upper and bottom layers are coupled since no assumption is 
made about the hypolimnion velocities -being negligible. 

The solutions are very dependent on the wind conditions. For a uniform wind stress 
of moderate strength, the "zeroth order" horizontal pressure gradients in the hypo- 
limnion are zero and the only hypolimnion velocities are in a thin boundary layer adja- 
cent to the thermocline. For this case, the bottom topography and value of the hypo- 
limnion eddy diffusivity only affect the solution to a small degree. However, for a wind 
stress with unit order gradient, small hypolimnion horizontal pressure gradients do 
occur producing significant hypolimnion geostrophic (inviscid) velocities. For this case, 
the value of the hypolimnion eddy diffusivity has a large effect on the thermocline shape. 
In addition, the shape of the lake bottom also becomes important. With unit order wind 
stress gradients, large differences in thermocline shape and horizontal velocities occur 
between the asymptotic solution which assumes small hypolimnion eddy diffusivity and 
previous solutions which uncouple the two layers by assuming the hypolimnion velocities 
are zero (hypolimnion eddy diffusivity very large). 


* Professor of Geophysics and Engineering, Case Western Reserve University, 
Cleveland, Ohio. Professor Lick’s work was supported by the National Science Founda- 
tion. / 



INTRODUCTION 


Observations of stratified lakes suggest that during steady wind periods they may be 
modeled by considering the lake to be made up of two homogeneous layers each with dif- 
ferent densities and eddy diffusivities with the interface between the two layers being 
located at the thermocline. This two-layer model is analyzed here. A more complete 
discussion of the applicability of the two- layer model is given in reference 1. 

Welander (ref. 2) and many others have studied a two- layer model for the oceans 
which included the Ekman dynamics assumption that the shear stress at the thermocline 
and the ocean bottom are proportional to the geostrophic velocities in either the hypo- 
limnion or epilimnion. Others such as Hamblin (ref. 3) have used this approach in the 
Great Lakes. In fresh water lakes, the thickness of the friction layer generated by the 
wind is of the order of or greater than the average thickness of the epilimnion; and there 
is no real geostrophic velocity in the epilimnion to which the shear stress at the thermo- 
cline can be made proportional. Therefore, Ekman dynamics cannot be used. In this 
analysis, an extension of the shallow lake equations originally derived by Welander 
(ref. 4) and shown by Gedney and Lick (ref. 5) to yield good quantitative results for T.aVp 
Erie during uniform temperature conditions is used. 

In their two-layer analyses, Welander (ref. 2) and Hamblin (ref. 3) initially un- 
coupled the bottom layer from the top layer by assuming that the motion in the bottom 
layer can be neglected (quasi-compensation assumption). After finding the position of 
the thermocline with this assumption they then allow small motions in the bottom layer. 
This small motion is calculated using the thermocline position and shear stress calcu- 
lated from the upper layer solution. 

As will be shown later in this report, the quasi- compensation assumption does not 
appear to be valid for a two-layer lake subjected to a spatially variable wind stress. 

We therefore make no assumption pertaining to the hypolimnion velocity magnitude. In 
order to solve the governing equations we assume that the ratio of hypolimnion to epi- 
limnion eddy diffusivity ^ v yi2^ v M\) * s sma ^ anc * perform an asymptotic expansion in 
that parameter. The solution is valid for any value of the vertical Ekman number (E v j) 
in the epilimnion. For the normal case in a large lake, E yl in the epilimnion is of 
order 1. For moderate winds, the value of v M2 which makes « 1 is only 

6 to 12 times that of the molecular kinematic viscosity. According to Hutchinson (ref. 6) 
and Sundaran, Esterbrook, Piech, and Rudinger (ref. 7) this is the correct magnitude 
f° r v y[ 2 ‘ Larger hypolimnion eddy diffusivities than those considered here have been 
investigated numerically by Gedney, Lick, and Molls (ref. 1) and Lee and Liggett 
(ref. 8). 

In order to determine the error induced by making the quasi- compensation assump- 
tion, the results of the shallow lake equations with the quasi- compensation assumption 
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included are given and compared with the asymptotic results of the complete two- layer 
equations. 


SYMBOLS 

Aj coefficient in eq. (7) 

Ag coefficient defined by eq. (8) 

a^ 

a 2 + 

B, C coefficients in equations (32) and (33) 

T >2 coefficients in eqs. (9) and (10) 

9D 2 /3£ 

dj epilimnion friction depth 

d 2 hypolimnion friction depth 

o 

E yl epilimnion Ekman number, "Ml/ f c L 

f arbitrary function 

f c Coriolis parameter 

f j unit order depth function 

f 2 unit order wind stress curl function 
G coefficient in equations (32) and (33) 

g acceleration of gravity 

H coefficient in equations (32) and (33) 

h nondimensional lake depth 

h m 9h / 9m 

h t 3h/0t 

K coefficient in eqs. (7) and (8) 

L reference length of lake 

Mti + iM yl , epilimnion volume transport 

M T2 ^x 2 + iM y2’ e Pili mnion volume transport 
m coordinate normal to boundary or curve 



m X ’ m y x " and y-components of m 
m unit normal vector to boundary 


d_ 

3n 

3 

3n* 


P 

s 



u,v,w 


X, y,z 



V 

P M1 

V M2 

k 

P 

P 3L 

P T 

Ap 


— + i — 
3x 3y 



3x 3y 

nondimensional pressure 
integer exponent 

coordinate tangential to boundary or curve 
wind velocity 

dimensionless velocity, respectively, in x-, y-, z-directions 
dimensionless Cartesian coordinates 

2 '«rsf/ d l 

2»5 re f/ d 2 

U 1 + ^ V 1 
u 2 + iv 2 
V<Vr> 

boundary layer stretching parameter 
nondimensional surface displacement 
? + £ 

epilimnion eddy diffusivity of momentum 
hypolimnion eddy diffusivity of momentum 
nondimensional thermocline depth 
density of water 
density of air 
pi/p 2 

<p 2 - p 1 )/p 1 

r x + ir y ’ dimensionless wind shear stress (r w /r^ ef ) 

(1 + i)/2 
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Subscripts: 

b value at lake boundary 

c value at lake center 

ref dimensional reference quantity 

x x- component 

y y- component 

0 value at x = 0, y = 0 

1 epilimnion 

2 hypolimnion 
Superscripts: 

(i) imaginary part 

(r) real part 

(0), (1), (2) zeroth, first, second term in asymptotic expansion 
— dimensional component 

~ unit order function 

+ stretched coordinate 

BASIC EQUATIONS AND BOUNDARY CONDITIONS 

In the present analysis, the lake is considered to be composed of two layers of dif- 
ferent density as shown in figure 1. In each layer, the basic approximations are that 
the water density is constant, the vertical eddy viscosity is independent of position, the 
pressure is hydrostatic, and the lateral friction and nonlinear acceleration terms can 
be neglected. Effects due to the Earth's curvature and to the variation in Coriolis force 
with position are neglected since the length scales of lakes are much less than the radius 
of the Earth. The explanation for the density and eddy diffusivity being considered con- 
stant in each layer is given in reference 1. Gedney and Lick (ref. 5) have shown the 
other assumptions to be good approximations for the Great Lakes. 

The analysis performed here will assume a steady wind. As is well known, the 
location of the thermocline in a stratified lake being acted on by a steady wind will 
slowly become deeper . The rate of deepening has been measured in Cayuga Lake (see 
ref. 7) to be in the neighborhood of 15 to 30 centimeters per day. This rate is of such 
magnitude that the time derivative inertia terms in the momentum equations can be 
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neglected and the problem is then "quasi- steady. " With the steady wind restriction 
plus the other assumptions stated, the nondimensional continuity and momentum equa- 
tions applicable to each layer shown in figure 1 are as shown in reference 1: 

3u 1 3v, 3w, 

— - + — - + = 0 (1) 

3x 3y 3z 


3u 0 3v 0 3w„ 

— + — + — _ = 0 
3x 3y 3z 

1 a2r i i r _ 1 _ 1 ap l 

2 a_2 2 * 2n 3n 2n 3n 

a ^ 3z 

2 

J_ 5 r 2 _j. r _ p r /34 [ 3g\_ p r 9p 2 

2 . 2 2 2 27T \3n 3n/ 2n 3n 

0?2 dz 


( 2 ) 

(3) 


(4) 


The Cartesian coordinate system used (see fig. 1) has x increasing eastward, y north- 
ward, and z vertically upwards with the corresponding velocities being m, Vj, and Wj 
where j = 1 indicates the epilimnion and j = 2 indicates the hypolimnion. In these 
equations T ^ = u^ + ivj, r 2 = u 2 + iv 2 , £ is the lake surface, £ the thermocline posi- 
tion, p. (j - 1, 2) the pressure, 3/3n = d/dx + i(3/3y), o-j - 2 ^ ref / d 1 > = 27r £ re f/ d 2> 

and p r - p^/p 2 where dj and Pj (j = 1, 2) are the friction depth and water density in 
either the epilimnion or hypolimnion. In forming these nondimensional equations, we 
have used the following relations: 


L ^ref 

y=I U=JL 

L u ref 


^ref u ref 

«=-I- w= _^_ 

^ref w ref 
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where 


Lr 


w 


ref 


5 ref 


PjdiS 




^ref 


fref 

Ap 



TIT. 


W 


u 


ref 


ref 


p l d l £ c 


Ap =■ 


w. 


; ref 


ref 


U. 


ref 


and the overbar indicates the dimensional quantity. Here is the reference wind 

stress, g the acceleration due to gravity, f c the Coriolis parameter, (j = 1, 2) 
the vertical momentum eddy diffusivity, and L the reference length of the basin. 

The system of equations (1) to (4) must be solved subject to the following boundary 
conditions: 


_w _w . • _w _ 1 r ] 

T = T +1T = — 

x y a, 3 z 




at z = 0 


r 2 = w 2 = 0 at z = -h(x, y) 


r 1 = r 2 




> 


i 2 

1 d 2 1 


ar i do ■« sr , 


> at z=5(x, y) 


3z d 2 p r 3z 


( 5 ) 
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( 6 ) 


Wj = 0 at z = 0 




w 


= u 9i + v ii' 
3x 1 3y 


> 


> at z = 4 (x, y) 


w 2 =u 2-| + v 2f 

3x 3y> 


where h is the depth of the lake and r™ and rj are, respectively, the x- and y- 
components of wind shear stress. The boundary conditions on w 1 and w 2 at 
z = £ (x, y) are such that there is no flow normal to the thermocline. For strong enough 
winds the thermocline may intersect the surface of the water. This case is not consid- 
ered here. The details associated with the derivation of equations (1) to (6) can be found 
in Gedney, Lick, and Molls (ref. 1). 

We have said nothing about the average thickness of the epilimnion which is deter- 
mined by the thermal problem. In this report, the z value of the thermocline will be 
specified at one point on the boundary so that the average thickness of the epilimnion is 
consistent with observations. 


GOVERNING EQUATIONS FOR THE HORIZONTAL VELOCITIES 

Equations (3) and (4) subject to the boundary conditions (5) can be solved for the 
horizontal velocities in the hypolimnion and epilimnion as functions of the surface dis 
placement £, the thermocline position £, the lake depth h, and the wind stress r w . 
The resulting solutions are 


w «i((-z) wa.^+z) w wa.z . 

r 1 = A ie 1 +A ie 1 + — e 1 +LK 

w n 3n 


£ ^ z < 0 


( 7 ) 


r 2 . (a, - ^ i-A e tt ’“ 2<h+Z> - A 2 e'"“ 2<h+Z> ♦ 

\ n dnj 1 n d n 


-h < z < £ 


( 8 ) 


where 


A 1= ^i2(i + e 

1 jtK 3n \ 


-2a 9 -a 9 ^ 

* - 2e z 



■2a r 


a. 
w „ 1 

[1 + e + 

( / w K 


-2a r 


1 - e 


- 6 1 + e 


-2a r 
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A 2 ' 


w o a l" a 2 ip 


2t e x e 
co K 


r 

nK 


4 2a l . iVl - e"* 2 




M. . ±. e '% 2 * 1 

3n 7rK \ 


iif^ 

3n 


= coo!^ 

&2 = + 4) 


CO = ■ 


1 + i 


5 =■ 


d lP r 


in =i£ + ii 

3n 3n 3n 


K = fi/l + e~ 2a2 Yl ♦ - (l - e 2 % 231 - l) 


The Wj and w 2 velocities can be determined by differentiating equations (7) and (8), 
substituting in equations (1) and (2), and then integrating. 


ASYMPTOTIC EXPANSION OF TWO-LAYER LAKE EQUATIONS 

As explained in the INTRODUCTION we consider here the case when j/ M2 / ^Ml ^ 1 

so that 6 = (l/p r ) y v m 2^ v MI K< X ^r = p l //p 2 ~ lm °)‘ We then P erform an asymptotic 

expansion which is to be valid as 6 — 0. The asymptotic expansion which we perform is 
valid as 6 — 0 provided that two conditions are met. The first condition is that 

-a 2 -2ttcop (h+4)/5d 1 3 

e 1 = e r 1 < 6^ 


In the limit as 6 — 0 this requires that h + 4 * 0 or that the thermocline does not in- 
tersect the lake bottom (note in the analysis h is always positive and 4 is always neg- 

/ 2a. i \ // \ 

ative). The second condition is that e + 1) /je - 1] is of order unity (a^ = c va^), 
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which requires that the thermocline be below the surface of the lake by a certain amount 
(i.e. , | < 0). The bounds on h and £ which ensure that these two conditions are met 
for the cases to be calculated will be discussed subsequently. 

The cases to be calculated will be for - 

p 

(1) A moderate wind of 5.2 m/s = 16.8 cm /sec) 

(2) 1/576 == v -^/ < 1/144 (^ M2 is therefore a few times larger than the 

molecular diffusivity) 

(3) p 2 - p x = 0. 002203 g/cm 3 

(4) Lake length of 96 km 

(5) Wind shear stress r^ ef st l dyne/cm 2 

For these conditions the values of 6, a 2 , and are 

— < 6 < — 

24 12 


84 < oi 2 < 168 


<3^ =8.4 


~ a 2 

With this range of 5 and a 2 , e <6 
10 percent deeper than the value of | £ j. 

/2a i 

ters. For =8.4, the quotient e 


providing the lake depth is approximately 
The value of £ ref is approximately 30 me- 


+ 1 


2a! 


- 1 is of order unity provided the 


minimum values of |£ | are of the order of 3 to 4 meters. 

With these restrictions understood we now proceed with the asymptotic expansion. 
First, 



The expansions for and A 2 are 
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^922 + ii^ + 0(6 3 ) 

ji dn n 9iy 



ip r 9r? | i 9^ 
n 9n 7r 9n 


+ ^in + o(5 5 ) 

jr 9n 


When the expansion for A 2 is substituted into equation (8), some cancellation takes 
place to give the more compact form 


r 2=' 


9 w o a * 

-2t e 


9n 7T 9n- 




iH e " WQ! 2 (h+z) 
it 9n 


+ ^£^H + 0(5 2 ) (8a) 

7T 9n 


When Aj is substituted into equation (7), no such simplification for results. 

The velocity equations (7) and (8a) can be integrated vertically to give the volume 
transport equations 
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r Tl = M vl + iM v1 = f° «i r 1 dz = -2ir w _ iflii I 
y •'l it 9n 


- 6 1 + a 


J T2 " M x2 + iM y2 “ 


L p - 


id Z =^^ (h+{ )-^^ 

it 3n i 3n 


+ 31 + 5 


Integrating equations (1) and (2) in the z-direction and applying equation (6) give 


Real/ = 0 


where 


Real 



= 0 


( 12 ) 



and 

D 2 (|) = 


Terms of order 6 3 have been neglected in equations (9} and (10). The quantity M™. 
(j - 1, 2) is the total volume transport in either the epilimnion or hypolimnion, and 
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M - and M • are, respectively, the volume transport components in the x- and y- 
xj yj 

directions. 

Equations (9) to (12) can be combined to give the governing equations for £ and £. 
This can be most readily accomplished by forming 


Real 




- 0 


which results in 


6p r V 2 r) + 


/jty 0h _ 0 tj 0h\ 
\0y 0x 0x 0y/ 




+ 0 ( 6 3 ) 


(13) 


Combining equations (10) and (12) results in the second governing equation 



3 

Equations (13) and (14) which neglect terms of order 6° are to be solved for the depend- 
ent variables subject to the boundary condition that the volume transport normal to the 
boundary is zero. This boundary condition is then 


M xl m x + M yl m y = 0 


(15a) 


M x2 m x + M y2 m y = 0 


(15b) 


where m x and my 


are the x- and y- components of unit normal m to the boundary. 


SOLUTION OF GOVERNING EQUATIONS 

The solution of the governing equations (13) to (15) will be performed using an 
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asymptotic expansion procedure . To develop the asymptotic expansion, we substitute 


7j =r? (0) + 6r? (1) + 5 2 r/ 2) + . . . 

4 = + 5 2 ^ + . . . 

5 = + 6 2 ^ 2 ) + . . . 

p 2 - p^ + 5p^ + 6 2 p^ + . . . 




( 16 ) 


into equations (13) and (14) to give 

5p r vy°> + 4/,® + s a ,| 

\ 3y 3x 3x 3y / l \ 


3 rj^ _3h _ dr,W 3h\ 
3y 3x 3x 3y/ 


V 3y 3x 3x 3y/ \3x 3y / 


vV°) + v 2 v {1) + (V 0) 9h _ 3r^ih\ , a l 5 (dy 

ir n n y 3y 3x 3x 3y/ 77 \ 3 


^ 3h _ 3-r/ 1 ^ 3h^ 


3y 3x 3x 3y/ 


a i /s„< 0) 3f (0) 3^t°) es (°>\ , aT) d) 3 g(0)-\ 

7T \ 3y 3x 3x 3y / 77 \ 3y 3x 3x 3y / 


+ ^< V°> 3g (1 > . si (0) a^ (1) V 8Eeal J a 


77 \ 3y 3x 3x 3y 


3n* 


t W D 2 (|<°)) 


„ Sr)®* , 2u S?*°* 

Pj, + 

77 3n 77 3n 


2a> 


1 -5 2 Reali-^~ 
3n* 


-v 


2u> 

— Pr 


3r? 


(1) 


3n 


2cu 3? 


(1) 


77 3n 


+ 0(6^)(terms involving first derivatives) 


+ 0 ( 6 ' 3 ) = 0 


(14a) 


14 



Similarly, the expansion form of the boundary conditions can be obtained by combining 
equations (9), (10), (15), and (16). 


Case 1 - Uniform Wind and Constant Depth 


Consider the asymptotic solution for a constant depth square basin which is acted 
on by a uniform wind over the entire water surface. As determined. from equation (13a), 
the zeroth order term (the term which is of order 6) from differential equation (13) is 

vV°) = 0 (17a) 


The boundary conditions for are determined from the zeroth order terms of equa- 
tion (15b) and are 


, (0) 

^22 — = o at x = 0, 1 
3y 

► 

a (0) 

^22 — = o at y = 0, 1 
3x 


(17b) 


and, therefore, t/ 0 ^ = constant on the entire boundary. 

As is well known, the solution of the boundary value problem represented by equa- 
tions (17) is = constant. Therefore, to zeroth order, the horizontal pressure gra- 
dient (see eq. (4)), 


a P2°) _ 3t/°) _ 3^°) [ 3g(°) 
3n 3n 3n 3n 


in the hypolimnion is zero. This does not mean the velocities in the entire hypolimnion 
are zero. There is a shear stress transfer at the thermocline which imparts motion to 
the hypolimnion water in a thin boundary layer adjacent to the thermocline. Below this 
boundary layer, the zeroth order velocities in the hypolimnion are zero. This can be 
shown from equation (8a) using a typical value for of the order of 168. 

It is possible to obtain £ ^ along the boundary from the zeroth order component of 
the equation (15a) boundary conditions. This results in 


/ 15 



(18a) 


2T W -^iil^5<°>=0 

' v 7 t dy 

2t w - — = o 

IT dX 


■N 

at x = 0, 1 

> 

at y = 0, 1 


where we have used the relations 


34 (°) _ _ 3 g(°) 

3x 3x 

and 

34 (°) _ _ 3 ^°) 

3y 9y 

These equations can be directly integrated for the value of 4 ^ along the boundary to 
give 






at x-0 


4 


( 0 )_ 




at y = 0 


>- 


4 


(°) = 



at x = 1 


(18b) 


4 


(°) = 



at y = 1 


j 
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where is the value of £ ^ at x = 0 and y = 0. 

Since is a constant, the zeroth order component of equation (14a) yields nothing. 
The second-order terms of equation (14a) give a governing equation for £ but it in- 
volves Therefore, to calculate £ we must also solve for The governing 

equation for ij^ is obtained from the first-order terms of equation (13a). These 
and £ ^ governing equations are 


vV 1 ) = 0 


(19) 


v 2 | <°> + tt fo (1) ag (0) _ja^ + J_ t w 

\ ay 3x ax 3y J 


x 




dx 


9y 


_w 


+ T.. 


D&(°))ii 

^ dx 


( 0 ) 


D (r)(^(°))ii 


( 0)1 


9y 


= 0 (20) 


where D 2 ^ = 3D 2 /2£ and the superscripts r and i indicate, respectively, the real 
and imaginary parts. In deriving equation (20), we have used the fact that 


3t 7 (0) _ 3£ (0) , ag (0) = 0 

3n 3n 3n 

These coupled equations are solved using the £^ boundary conditions specified by 
equation (18) and the following boundary conditions derived from the first-order 
terms from equation (15b): 


Real 


r w D 2 (^°))] + l(ii 


(0) g^(9)\ a -\ anU) 


1 dr}2 


9y 


dx 


9y 


(h+ £ (0) ) = 0 at x = 0,1 


(21a) 


jW[r w D 2 (£^) 


m (h + £ (°)) = 0 at y = 0, 1 (21b) 

3x 3y ) n dx 


The solution for £^ and were obtained numerically using finite differences. The 
numerical method used was Newton's method with point under relaxation. 


17 



Case 2 - Uniform Wind With Variable Depth 


As has been just discussed for case 1, a uniform wind acting on the water surface 
of a constant depth basin produces no zeroth order flow in the hypolimnion except for a 
thin boundary layer adjacent to the thermocline. Implicit in the assumptions of this 
whole analysis is that this thin hypolimnion boundary layer does not contact the lake bot- 
tom. As a result, there are no flow velocities in the region of the lake bottom which 
could be affected by variations in depth. Therefore, variations in bottom depth should 
not appreciably effect the zeroth order solution for the uniform wind case. This conclu- 
sion has been verified by reference 1 where the complete two-layer lake equations have 
been solved numerically for a lake whose depth is constant along the shore but is vari- 
able in the interior . 

In this section the asymptotic solution rj ^ for a variable depth lake will be obtained 
to again demonstrate this conclusion. Here we consider a lake whose depth is described 
by 


h = h fe + 16. 0(h c - h b )(x - x 2 )(y -y 2 ) 0 < x < 1, 0<y<l (22) 

where h^ is the depth at the lake boundary which is taken to be a constant and h c is 
the maximum lake depth which is taken to occur at the center of the lake. Note that the 
h derivatives tangential to the boundary are directly proportional to the distance from 
the boundary. At a small distance A from the boundary, the tangential depth deriva- 
tives are of order A. 

The zeroth order terms from the governing equation (13a) is 

(23a ) 

3y 3x dx dy 

where we assume at this point in the analysis that no boundary layer exists which 
would make the V 2 tj^ term in equation (13a) or order one. The general solution of 
equation (23a) is 


t/°) = f[h(x,yj] (23b) 

where f is an arbitrary function of the depth h(x, y). The boundary conditions on 
can be shown from the expansion of equation (15b) to be = constant. Since by equa- 
tion (22) h is constant on the boundary, equation (23b) satisfies the boundary condition 
for „(°). 

In order to determine the form of f we investigate the zeroth order term of equa- 
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tion (14a) which is 


9t/°) 0h dr/^_ 5h [ dy^ d^ dy^ 3|j^ = Q 

0y 0x 0x 0 y 0y 0x dx dy 


(23 c) 


Substituting equation (23b) into this equation results in 

0h0|^_0h0|^ =o (23d) 

0y 0x dx dy 

which has the general solution 

£ (0) =f[h(x,y)] (23e) 

where f is again some arbitrary function. The boundary conditions for £ are deter- 
mined from the zeroth order term of (15a) and are the same as equation (18b). Since h 
is constant on the boundary, no function f will be able to satisfy these boundary condi- 
tions. 

A possibility exists that there is a boundary layer so that equation (23d) does 
not apply in the boundary region. To determine the boundary layer equations, it is nec- 
essary to introduce the stretched coordinates 


m + = — , t + = y along x = 0, 1 
e(6) 

> 



along y = 0, lj 


(24) 


in equation (14a). The resulting zeroth order boundary layer equations are then 


5 0V 0 ) 

e 2 0m +2 


+ O', 


3f 

0h 


(-1) S+1 0g (O) 

_ e 0m + 


ht + (-l)’ 


S 0|^ h 


m 


= 0 


d€ 


(25a) 


where this equation applies along x = 0, 1 when S = 1 and along y = 0, 1 when S = 2 
and where h and h. are, respectively, the normal and tangential partial derivatives 
of h. All derivative terms except for h^. and d% K ’/d t in this equation have been 
scaled such that they are of order one. Equation (22) shows that in the boundary layer 
region 
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(25b) 


h t =^- = €f 1 (m + ,t + ) 
3t + 


where f^ is of order one along the entire boundary. In order to estimate the magnitude 
of 0^°V3t + in the boundary layer we must consider the boundary conditions (eq. 18(b)) 
and the interior solution (eq. (23(e)). The boundary conditions show that 3£ ^7 3t + will 
be of order one along y = 0, 1 but zero along x = 0, 1. The interior solution (eq. (23(e)) 
shows that near the edge of the possible boundary layer 


31; _ df 3h . df * ,_+ *+ 


3t 


dh 


3t 


= € ~ fj(m ,t ) 
+ dh 1 


(25c) 


where e is the thickness of the boundary layer. A Taylor series expansion for £ ^ 
using equations (18b) and (25c) will show that d£^/dt + in the x = 0, 1 boundary layers 


will be of 0(e). In the y = 0, 1 boundary layers 3£ ^/dt + is of 0(1). 

When equation (25b) is substituted into equation (25a) and the estimate is used for 
94 ^/dt + , the boundary layer equation along x = 0, 1 becomes 


6 (0) gt(0) , , 

--2-5 — + a* — f, (m + ,t + ) + O(e) - 0 

£ 2 8m +2 3h 8m + 


When the principle of least degeneracy is used, the proper scale is chosen as e = 6^^ 
and the final form of the x = 0, 1 boundary layer equation became 


3V°> 


3m 


+2 


+ a. 


3f 3£ 


( 0 ) 


f 1 (m + ,t + ) + 0(6 i/A ) = 0 


9h 3m + 


l/2\ 


The boundary layer resulting from this equation is an increasing or decreasing exponen- 
tial depending on the sign of 3h/3t + = ef 1 (m + ,t + ). Only exponential decreasing functions 
can be allowed. Because 3h/3t + = ef^ changes sign about y = 0. 5, it will only be pos- 
sible to have a boundary layer along one -half of the x = 0, 1 boundaries. With equa- 
tions (22) and (23e), the values of outside the boundary layer will be symmetric 
about y = 0. 5. But this means that a % ^ boundary layer will be needed along the en- 
tire x boundary. We conclude from this that no boundary layer is possible along 
the x = 0, 1 boundaries which will permit the £ boundary conditions to be satisfied. 
As a result, no solution exists when we assume a nonconstant form for equa- 
tion (23b). The only other possibility is that 77^ is a constant. 
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Before finally concluding that 77 ^ is constant for this variable depth basin^case, 
we must admit the possibility of an 7 ] ^ boundary layer that would make the V 77 
term in equation (13a) of order one. If this is done, the same arguments just used to 
show that no boundary layer existed will show that no 77 ^ boundary layer can 
exist. Therefore, we conclude that 77 ^ is a constant and that the zeroth order horizon- 
tal pressure gradient in the hypolimnion ^3p^ V 3n - dr] ^ / on) is zero whether or not the 
lake depth is variable. 


Case 3 - Wind Stress Gradients of Order 6 With Constant Depth 


The zeroth order terms from the governing equation (13a) when the lake depth is 
constant but the wind shear stress curl (arj'/sx - 3T*/3y) is of order 6 is 

vV 0 >=2£ f (26) 

p r \zx ay/ 


where 



3 w\ 

8r y _ dr x 

3y 


is of order one. The boundary conditions for 77 ^ can be derived from equation (15b) 
and again 77 ^ is a constant on the boundary. This boundary value problem for 77 ' 
has a nonconstant solution which will result in hypolimnion pressure gradients 
(dp^/dn = 377 (° Van ). 

The governing differential equation for £ w is obtained from the zeroth order com- 
ponent of equation (14a). This equation for is singular and on the interior region of 
the basin is 


3^(0) ^(Q) ^<0) __ p 

3y 3x dx dy 


(27a) 


In order to satisfy the boundary conditions, a "boundary layer" must exist along por- 
tions of the boundaries. To determine the boundary layer equations, it is necessary as 
in the previous case to introduce stretched coordinates given by equations (24). The 
resulting zeroth order boundary layer equations are then 
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6 syp) 

e 2 3m +2 


+ a , 


(-irl t| WM^ +( . 1)S l,(0)3 i 


( 0 )' 


3t 


3m + J 


= 0 


where this equation applies along y = 0, 1 when s = 1 and along x = 0, 1 when s = 2 
and where and are, respectively, the normal and tangential partial derivatives 
of 7} ^ given by the solution of equation (26). Because of the boundary condition on t/°) 
(i. e. , 7j 0 = constant on the boundary), tj^ within a possible boundary layer thickness 

e can be^shown by a Taylor series expansion to be O(e). We denote here = «4 0) 

where 77^ is of order one. Using the principle of least degeneracy the proper value of 
e is then e = 6^ and the boundary layer equation becomes 


3V°> 


3m 


+2 


+ ai-, 




(or 


3t 


3m + J 


= 0 


(27b) 


Instead of obtaining the interior and boundary layer solutions separately, we obtain the 
complete 77^ solution from the combined equation 

6 vV°> + a . ia® 0 (2 

\ 3y 3x 3x 3y J 

The boundary conditions for % are the same as equation (18). The solution of equa- 
tion (28) is easily obtained numerically using successive overrelaxation. 


Case 4 - Wind Stress Gradients With Variable Lake Depth 

The results to be given subsequently for case 3 where the basin depth is constant 
and the wind stress curl is of order 6 show that significant geostrophic velocities occur 
in the hypolimnion due to the finite values of 3p^ V 3n = drj^/dn. With significant hypo- 
limnion velocities, variations in lake depth should be very important. Consider again a 
square lake whose depth is described by equation (22) which is acted on by a wind with a 
stress curl (W^/Sy - drj/dxj of order 6 (case 4a) and order 1 (case 4b). At the lake 
boundary, normal derivatives of the wind stress are nonzero but tangential derivatives 
are assumed zero. 

Case 4a - wind stress curl of order 5 . - The zeroth order governing equation for 

7} and for this case are the same as case 2. The solution is therefore rj^ = 
constant. The variation in lake depth is then enough to eliminate the hypolimnetic 
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horizontal pressure gradient (3i/°Van) which was created by the 6 wind stress curl 
when the depth was constant. The effect of variation of lake depth is therefore very 
strong. 

Case 4b - wind stress curl of order one . - For this case the zeroth order term of 
the governing equation (13a) is 


drj^dh _ a^Sh = 2rr_ f^y_ . _^x\ (29 

ay ax ax ay y 3x ay J 

The boundary condition on from equation (15a) is i/ ^ = constant. At the bounda- 
ries, equation (22) shows that the tangential depth derivative \ is zero and the normal 
depth derivative h is finite. Using this fact in the nonhomogene ous equation (29) 
shows that the tangential t/ 0 ^ derivative at the lake boundary is nonzero. As a result, 
equation (29) can not satisfy the = constant boundary conditions. We must there- 
fore consider an r/ 0 ^ boundary layer which will be compatible with a solution of equa- 
tion (29) in the interior of the lake. 

To determine the boundary layer we introduce stretched coordinates given by equa- 
tion (24) into equation (13a) which gives 



aV°> 


3m 


+2 


+ a * 


/ i\S+l p„(0) s 

til i3 — h t + (-l) s ia — h m 


3m 


at + 



where this equation applies along y = 0, 1 when s = 1 and along x - 0, 1 when s 2 
and where h t and h m are, respectively, the normal and tangential partial derivatives 
of h with respect to the boundary. Now equation (22) shows that 

h t = ef 1 (m + ,t + ) 

where i 1 is of order one. Again using the principle of least degeneracy gives the prop- 
er scaling as e = 6 1 / 2 . The boundary layer equation is then 


, iV 0> 

r 8m +2 


+ as 


(-l) s+1 fj(m + ,t + ) + (-l) s ^ h m (m + ,t + ) 


3m 


ar 



Instead of obtaining the interior and boundary layer solutions separately, we choose to 
solve the combined equation 
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(30) 


6 P r vV 0) + 


w 


a * 


3h _ dT]^ 3h\ _ 2n / dr y _ 
\ 3y 3x 3x 3y/ \ 3x 


,w\ 


dr 


3y 


over the entire region. The boundary conditions on 77^ are shown from the expansion 
of equation (15b) to be 77^ = constant. 

The solution of equation (29a) will result in 77 ^ - r/°\x,y) with a boundary layer of 
order 5 ' . The governing equation for 4 ^ is obtained from equation (14a), and the 
boundary value problem is easily seen to be singular. In the interior region, 

/ W 0) 3h 3 t? (Q) 3h \ + /3 t 7 (0) 34 (0) ag(0) \ Q 

\ 3y 3x 3x 3y/ \ 8y 9x 3x 3y / 

The equation valid in the boundary layer can also be readily derived. Introducing the 
equation (24) stretched coordinates into equation (14a) results in 


( 2 P r - l)n 


( 0 ) 

mm 



+ a 


(-1) S+ V°) 9€ (0) + (-1) S _(0) 3^ 

v v n m + T/ t — 


(0) 


3t 


3m 


. « a 2 t (0) , 


£ 2 9m +2 


where Tjj^ is the second partial derivative of t/ 0) with respect to m and is of order 
one or less. From the rj' 0 ' solution determined from equation (29a), 77^ (within a 

distance 6 1 / 2 from the boundary) was found to be where 77^°^ is of order 

one (this same behavior could also have been found by expanding 77^ in a Taylor series 
about the boundary). Two possible 4^ boundary layers can exist depending on the 
magnitude of 77 m in the boundary region. In some boundary regions an 77^) boundary 

layers exists so that 77^ = 77^/6 1/2 ; in other boundary regions v (0) = v (0) ( v (0) is of 

order 1). Where an V boundary layer exists, the 4 ^ boundary layer is of order 
6 (e = 6 / ) and the appropriate boundary layer equation is 


^“i<-i> s+1 >4 0,i ^=o 
+2 1 m 


3m 


3t 


In regions where no 77^ boundary layer exists, the 4^ boundary layer 
5 (e = 5 / ) and the appropriate boundary layer equation is 


is of order 
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^7 * to r ' + “l 

3m 

,( 0 ) 




at 


Sm 


= 0 


However, as with the solution for t; , we do not choose to obtain the interior solution 
and boundary layer solution separately but to solve the combined equation 

5 vV°> + 6 ( 2p - i)vy°> + a ^ jhV a (^ 1 . m < 0) $ 3 , o 


3y 3x 3x 3y, 


Again the boundary conditions on 4^ are the same as for equation (18). 


3y 3x 3x 3y 


(31) 


SOLUTION WITH ZERO VELOCITY IN THE HYPOLIMNION 


As discussed in the INTRODUCTION, the "quasi-compensation” assumption is com- 
monly made in two- layer lake analyses to find the thermocline position. The quasi- 
compensation assumption is that the horizontal pressure gradient and the velocity to be 
zero in the hypolimnion. Taking the hypolimnion velocity equal to zero is equivalent to 
taking the hypolimnion eddy diffusivity (y^) to be infinite. The zero hypolimnion veloc- 
ity case is considered here so as to be able to contrast its results with the results for 
the more reasonable case where the hypolimnion eddy diffusivity is small. 

The advantage of using the zero hypolimnion velocity assumption is that it results in 
a single governing equation for the thermocline position 4- If the boundary condition 
Wj = r i = 0 at z = 4 is used at the thermocline, the following governing equation and 
boundary conditions for 4 can readily be derived using equations (1), (3), (11), and 
(15a): 



= +2tt 





3B_3J_ + _3C _3Aw + /3C _3| _ 3B 3|\w 

34 3x 34 3y/ y \34 3x 34 dy/ y 

— 


(32) 


oii - H-2£ = -2 jt(btJ - CtT) at x = 0, f 

3x 3y v x y/ 

H ii + - -2 jt(CtJ + Brjf) at y = 0, 1 

3x 3y v y/ J 


(33) 
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where G, H, B, and C are given by 


sin ( oi-t £ ) - sinh (q'-.O 

G = + — 

cos (oi^) + cosh (a^) 


sin ( a, £ ) + sinh ( a* £ ) 

H 1 = — + a,* 

cos (a^i;) + cosh (a^, ) 


2 sin (a^/2) sinh (aj^/2) 
cos (q!j 4) + cosh (a^) 


2 cos {a*Z/2) cosh (a-,^/2) 

C = ± 1 

cos (a^i;) + cosh (q!j4) 

The derivation of these equations may also be obtained directly from the results given by 
Welander (ref. 4) for a shallow homogeneous lake. When h is replaced by -i- and 
3C/3x by -3£/3x, etc, Welander 's equations (11) and (13) lead to equations (32) and (33) 
after they are nondimensionalized in the manner given herein. It should be noted that 
equations (32) and (33) are not part of an asymptotic expansion but include terms of all 
order. Equations (32) and (33) can be solved numerically using successive underrelaxa- 
tion with Newton's method incorporated. At each iteration the equations are linearized 
about the £ values determined at the previous iteration. 


RESULTS 

The governing equations derived in the previous sections were solved for the posi- 
tion of the thermocline the position of the lake surface £, and the hypolimnion and 
epilimnion horizontal velocities u and v. All results are presented for a square basin 


with the following parameter values: 

Basin length, L, km 96. 6 

Epilimnion temperature, °C 22 

Hypolimnion temperature, °C 4 

Density difference, A p 0. 002203 

Density ratio, p 0.9977969 

*■ 9 

Epilimnion eddy viscosity, i > M1 , cm /sec 16.8 

Nominal wind velocity, U , m/sec 5.2 
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The value corresponds to that generated by a wind velocity of 5. 2 meters per 

second. Small variations in wind velocity about the nominal value will be allowed. The 
wind shear stress will be determined by 

tJ - -(0.00273)(p a )(U w ) 2 



where p a is the density of the air and the 0. 00273 wind drag coefficient is that proposed 
by Wilson (ref. 9). 

The wind drag coefficient and values used here are the same as the ones used 

by Gedney and Lick (ref. 5) for current calculations in Lake Erie during the uniform 
water temperature period (winter period). With these values agreement between the 
Lake Erie calculations and measurements were very satisfactory. The epilimnion in 
our two- layer model is considered to be of uniform temperature and for the cases cal- 
culated here will have an average thickness of approximately 20 meters. Since 20 me- 
ters is also near the average depth of Lake Erie, the value should be very similar 

to that used in the Lake Erie calculations. 

In obtaining all the thermocline solution the boundary conditions (18b) were adjusted 
to give the desired mean thermocline depth. 


Case 1 - 5. 2 Meter per Second Uniform Wind and Constant Depth 

Zeroth order asymptotic solution . - The zeroth order solution is governed by equa- 
tions (17) to (21) and gives the pressure gradient in the hypolimnion as zero; that is, 

dP2^ a^(0) ^ 3g(0) _ Q 
9n 3n 3n 

Therefore, there is no zeroth order geostrophic inviscid current in the hypolimnion. 
There are, however, hypolimnion velocities created by the shear stress at the thermo- 
cline. The zeroth order thermocline (£^) contours are shown in figure 2(a). The 
thermocline is at a shallow depth (4.6 m) at the upwind end of the lake and a large depth 
(=30 m) at the downwind end. There is very little tilting of the thermocline in the cross 
wind direction and, as determined by the equation (18b), there is no tilting at the upwind 
and downwind boundaries of the lake. The dimensional form of equation (18b) for the 
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thermocline boundary values is 


7 ( 0 ) _ 




g(P 2 - P x ) 


t(0) 

^0 


at y = 0, 1 


7 ( 0 ) _ 


r(o) 

■0 


at x = 0 


7 ( 0 ) _ 


g(P 2 - Pj) 


7 ( 0 ) 


at x = 1 


where the overbar indicates dimensional quantities. It is important to note that the 
boundary values of the zeroth order solution for the thermocline do not depend on (1) the 
eddy diffusivity values of (2) the depth of the lake. Furthermore, the length to width 
ratio of the lake has no influence on the £ ^ boundary values. This zeroth order solu- 
tion, of course, is only valid as long as dg/dj = ^ M2 /^ M1 « 1 and the lake is deep- 
er than the thermocline by the amount d 2 , which is the friction depth in the hypolimnion. 

In figures 2, the horizontal velocities are given at the lake surface and at depths of 
6.7, 15, and 20 meters. These velocity results have been determined from equations (7) 
and (8) assuming d 2 = 1. 5 meters. The dashed line included on some of the velocity 
plots is the intersection of the thermocline with the horizontal plane at the particular 
depth. The velocity patterns down to 10 meters are similar to those which occur in a 
homogeneous lake. The Coriolis force causes the deflection of the surface velocities to 
the right of the wind and a clockwise rotation of the current vector with depth. The ef- 
fect of the epilimnion thickness being smaller at the upwind end of the lake is easily 
discerned from the figures. At depths closer to the thermocline, the epilimnion veloc- 
ities rotate even more clockwise until they are in a southerly cross wind direction. As 
we go deeper to a region below the thermocline the velocities become very small in mag- 
nitude. At a distance approximately equal to d 2 below the thermocline the velocities 
are essentially zero. 

Effect of increasing the hypolimnion eddy diffusivity . - The effect of increasing the 
eddy diffusivity on the thermocline position is shown by considering solutions for values 

of 5^= 1/ p r (VW^)) °* 0* 0- 5 and °° for the uniform wind and constant depth 

condition. The 6=0 case is obtained from the zeroth order solution governed by equa- 
tions (19) and (21). The results for the 6 = 0. 5 case are obtained from the n um erical 
solution of the complete two-layer lake equations by Gedney, Lick, and Molls (ref. 1). 
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The results for the 6 = 00 case are obtained from equations (32) and (33) which assume 
zero hypolimnion velocity. As previously discussed in the INTRODUCTION, the zero 
hypolimnion velocity assumption has been used by others (refs. 2 and 3) to solve for the 
position of the thermocline since it uncouples the lower layer from the top layer. 

The thermocline positions for 6=0, 0.5, and °° are shown in figure 3. In each 
case the minim um depth of the thermocline is approximately 8. 4 meters and occurs 
along the right boundary. As 6 is increased (i/j^ increased) the thermocline tilt in- 
creases though not dramatically. For uniform wind conditions, the 6 = 00 case, which 
assumes zero hypolimnion velocities, does not differ greatly from the more realistic 
6 < 1 case. 


Case 2 - Uniform Wind With Variable Depth 

As has been shown by the case 1 asymptotic results with 6 « 1, a uniform wind 
acting on the water surface of a constant depth basin produces no zeroth order flow in 
the hypolimnion except for a thin boundary layer adjacent to the thermocline. The 
zeroth order horizontal pressure gradient for case 1 is zero. From these results we 
expect that a variation in lake depth should produce no or a negligible effect when the 
winds are uniform. This was shown in the analysis section to be true. The zeroth 
order horizontal pressure gradient in the hypolimnion is zero for a uniform wind with 
variable depth. 


Case 3 - Wind Stress Gradients of Order 6 


In this section we examine the effect of a variable wind stress of the form 


rj = -0.914 - 8.0 6 (x - x 2 )(y - y 2 ) 


r™ = 0. 0 


0 s x < T 
0 ^ y ^ 1 


(34) 


acting on the surface of a constant depth basin. This t™ has the largest wind stress at 
the center of the lake which often occurs in actual lakes. (Note that x = 1, y = 1 corre- 
sponds to x = 96. 6 km, y = 96. 6 km. ) This wind stress also creates a maximum wind 


stress curl 


dr^/dy 


of 26 at x = 0. 5 and y = 0. 0 and 1. 0. The governing equations 
for this case are equations (26) and (28). 

The thermocline depth contours are shown in figure 4(a) for 6 = 1/24. These re- 
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suits show that the thermocline shape is very sensitive to spacial variations in the wind 
stress. Whereas the thermocline generally has very little variation in the cross wind 
direction when the wind stress is constant, a wind stress curl of order 6 creates a 
large variation. 

The boundary layer of width 6 1//2 L predicted by equation (27) occurs along 
several portions of the boundary as shown in figure 4(a). The "outer” solution governed 
by equation (27a) extends throughout the rest of the lake other than these boundary layer 
regions. 

Shown in figure 4(b) are significant horizontal velocities below the thermocline at a 
depth of 32 meters. These velocities occur below the viscous boundary layer adjacent 
to the thermocline and therefore are geostrophic (inviscid) velocities created by the 
nonzero horizontal pressure gradient in the hypolimnion. The governing equation for 
the geostrophic velocity can be obtained from equation (4) by dropping the viscous terms 
and is 


9p 2 ^ _ 3 ? _ 34 + 

3n 3n 3n 3n 


-771 


i4 0) 


+ IV, 


(0)1 
2 J 


The geostrophic velocity occurs at a right angle to the direction of the horizontal pres- 
sure gradient Vp^ = As a result, the geostrophic velocities are tangent to the 

*7 0 contours shown in figure 4(c). The geostrophic flow in the region of the lake bottom 
is, of course, brought to rest in a velocity boundary layer at the lake bottom. 

With significant geostrophic flow in the hypolimnion, we expect a depth variation to 
be very important. As shown in the analysis section, the inclusion of a variable depth 
with an order 6 wind stress curl eliminates the hypolimnion horizontal pressure gra- 
dient (i.e. , Tj is constant) and therefore the hypolimnion geostrophic velocity. 

The influence of bottom topography is therefore very strong. We next consider the 
important case where the wind stress gradients are of order one over a basin of variable 
depth. 


Case 4 - Wind Stress Gradients of Order One 


In this section we consider the two- layer solution for a variable wind stress of the 
form 


= -0.914 - 2. 0(x - x 2 )(y - y 2 ) 


0 < x < f 
0 < y < 1 



( 35 ) 
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acting on the surface of a variable depth basin of the form 

Q = 2.0(1 - 2x)(y - y 2 ) 

0X 

= 2.0(1 - 2y)(x - x 2 ) 

3y J 

(Note that x = 1, y = 1 corresponds to x = 96. 6, y = 96. 6 km. ) The wind stress creates 
a maximum wind stress curl (dr^/dy) of -0. 5 at x = 0. 5 and y = 0. 0 and 1. 0. The 
asymptotic expansion solution to this case is governed by equations (30) and (31). 

The thermocline depth contours are shown in figure 5(a) for a 6 = 1/24. The £' 
solution contains boundary layers along y = 0, y = 96. 6, and x = 0 boundaries. These 
boundary layer regions are more clearly shown in figure 5(b) which is a plot of the ther- 
mocline depth along the x = 48.3 km (x = 0. 5) station. Figure 5(b) shows the thermo- 
cline to have a large shape variation in the cross wind direction. 

The horizontal velocities at the surface and at 6. 7, 20. 0, and 35. 0 meters are 
shown in figures 5(c) to (f). The figures show the effect of the shape of the thermocline 
on the velocities. The thermocline shape produces the clockwise gyre shown in fig- 
ure 5(e). Below the thermocline at a depth of 35. 0 meters geostrophic velocities of 
order 1. 0 to 5. 0 centimeters per second occur. As explained in the previous section 
(case 3), these geostrophic (inviscid) velocities are created by and are perpendicular to 
the horizontal pressure gradient Vp( 0) = The geostrophic velocities are there- 

fore tangent to the contours shown in figure 5(g). The high geostrophic velocities 
near the boundaries are created by the boundary layers shown in figure 5(g). 

In contrast to the uniform wind stress solution, the unit order wind stress gradient 
solution is very dependent on the value of the hypolimnion eddy diffusivity. Figure 5(b) 
shows the shape of the thermocline at x = 48.3 kilometers (x = 0. 5) for the same applied 
wind stress for the cases where 6 = 1/24, 1/12, and «*>. The 5 = « case corresponds 
to the solution governed by equations (32) and (33) where the velocities on the entire 
hypolimnion are assumed zero. The thermocline contour plots for the 6 =1/12 and 
5 = oo cases are shown, respectively, in figures 5(h) and (i). In all three cases the 
minimum thermocline depth of 6. 9 meters occurs along the right boundary. As the 
hypolimnion eddy diffusivity (i> M2 ) increases (5 increases), the cross wind variation in 
the thermocline depth can be seen to become much less. 


0 < x < 1 
02y<l 


CONCLUDING REMARKS 

The steady-state, wind-driven circulation has been calculated in a stratified lake 
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composed of two layers having uniform but unequal densities and eddy diffusivities. The 
solution is obtained from the governing equations assuming the ratio of hypolimnion 
(lower layer) to epilimnion (upper layer) eddy diffusivities (v^/v M1 ) is much less than 
one. A value of this ratio much less than one is generally recognized as being true for 
lakes during the late summer stratification period. 

The results indicate that the thermocline shape is strongly dependent on (1) spatial 
variation of the wind (on the order of 1 m/sec) over the lake, (2) the variation in lake 
depth, and (3) the magnitude of the hypolimnion eddy diffusivity. 

When the wind magnitude varies spatially over a lake, significant velocities occur 
in the entire hypolimnion. Previous two-layer lake models have calculated the thermo- 
cline position by making an interim assumption that the velocity in the entire hypolim- 
nion is zero (this is equivalent to assuming the hypolimnion eddy diffusivity is infinite). 
The thermocline shape calculated in this manner is found to differ by a large amount 
from the asymptotic solution which assumes a small hypolimnion eddy diffusivity. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, January 24, 1973, 

160-75. 
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(a) Thermocline depth contours, £' 0 ' (m). 

Figure 2. - Zeroth order two-layer lake solution with uniform wind and 
constant depth (6 - 1/241. 
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(e) Horizontal velocities at 20 meters. 
Figure 2. - Concluded. 
























































